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The scattering of random surface gravity waves by topography of small amplitude, and 
horizontal scales of the order of the wavelength, is investigated theoretically in the pres- 
ence of a an almost uniform irrotational current. This problem is relevant to ocean 
waves propagation on shallow continental shelves where tidal currents are often signifi- 
cant. Defining the small scale bottom amplitude normalized by the mean water depth, 
i] = h/H , a perturbation expansion of the wave action to order rj 2 yields an evolution 
equation for the wave action spectrum. Based on numerical calculations for sinusoidal 
bars, a mixed surface-bottom bispectrum, that arises at order 77, is unlikely to be sig- 
nificant in most oceanic conditions. Neglecting that term, the present theory yields a 
closed equation with a scattering source term that gives the rate of exchange of action 
between spectral wave components that have the same absolute frequency. This source 
term is proportional to the bottom elevation variance at the resonant wavenumbers, and 
thus represents a Bragg scattering approximation. With current, the source term for- 
mally combines a direct effect of the bottom topography with an indirect effect of the 
bottom through the modulation of the surface current and mean surface elevation. For 
Froude numbers of the order of 0.6 or less, the bottom topography effects dominate. For 
all Froude numbers, the reflection coefficients for the wave amplitudes that are inferred 
from the source term are asymptotically identical, as rj goes to zero, to previous theoreti- 
cal results for monochromatic waves propagating in one dimension over sinusoidal bars. In 
particular, the frequency of the waves that experience the maximum reflection is shifted 
by the current, as the surface wavenumber k changes for a fixed absolute frequency. Over 
sandy continental shelves, tidal currents are known to generate sandwaves with scales 
comparable to those of surface waves, with bottom elevation spectra that roll-off sharply 
at high wavenumbers. Application of the theory to such a real topography suggests that 
scattering mainly results in a broadening of the directional wave spectrum, i.e. forward 
scattering, while back-scattering is generally weaker. The current may strongly influence 
surface gravity wave scattering by selecting different bottom scales, with widely different 
spectral densities due the sharp bottom spectrum roll-off. 



1. Introduction 

Following the early observations of Heathershaw (1982), a considerable body of knowl- 
edge has been accumulated on the scattering of small amplitude surface gravity waves by 
periodic bottom topography. An asymptotic theory for small bottom amplitudes, that 
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reproduces the observed scattering of monochromatic waves over a few sinusoidal bars, 
was put forward by Mei (1985), leading to practical phase-resolving equations that may 
be used to model this phenomenon for more general bottom shapes (Kirby 1986). For 
sinusoidal bottoms of wavenumber I, Mei (1985) proposed an approximate analytical 
solution. In two dimensions (one horizontal and the vertical) this solution yields simple 
expressions for the wave amplitude reflection coefficient R, as a function of the mismatch 
between the wavenumber of the surface waves k and the resonant value 1/2, for which 
R is maximum due to Bragg resonance. Beyond a cut-off value of that mismatch, it was 
found that the incident and reflected wave amplitudes oscillate in space instead of de- 
creasing monotonically from the incident region. In three dimensions the Bragg resonance 
condition becomes k = I + k' and ui = ui', with lo and w' the wave radian frequencies 
corresponding to the wavenumber vectors k and k' through the linear dispersion relation. 

Other contributions have shown that higher-order theories are necessary to represent 
the sub-harmonic resonance observed over a bottom that is a superposition of two com- 
ponents of different wavelengths (Guazzelli, Rey & Belzons 1992). Such sub- harmonic 
resonance was found to have as large an effect as the lowest order resonance for bottom 
amplitudes of only 25% of the water depth, due to a general stronger reflection for rela- 
tively longer waves. However, these amplitude evolution equations are still prohibitively 
expensive for investigating the propagation of random waves over distances larger than 
about 100 wavelengths, and the details of the bottom are typically not available over 
large areas. Besides, a consistent phase-averaged wave action evolution equation is also 
necessary for the investigation of the long waves associated with short wave groups (Hara 
& Mei 1987). 

The large scale behaviour of the wave field may rather be represented by the evolution 
of the wave action spectrum assuming random phases. Such an approach was already 
proposed by Hasselmann (1966) and Elter & Molyneux (1972) for the calculation of wind- 
wave and tsunami propagation. A proper theory for the evolution of the wave spectrum 
can be obtained from a solvability condition, a method similar to that of Mei (1985) and 
Kirby (1988), but applied to the action spectral densities instead of the amplitudes of 
monochromatic waves. In the absence of currents the correct form of that equation was 
first obtained by Ardhuin & Herbers (2002, hereinafter referred to as AH) using a two 
scale approach. They decomposed the water depth H — h in a slowly varying depth H , 
that causes shoaling and refraction, and a rapidly varying perturbation h with zero mean, 
that causes scattering. This equation is formally similar to general transport equations for 
waves in random media (e.g. Ryzhik, Papanicolaou & Keller 1996), although the waves 
considered here propagate only in the two horizontal dimensions. The resulting scattering 
was shown to be consistent with the dramatic increase of the directional width of the 
wave spectra observed on the North Carolina continental shelf (Ardhuin et al. 2003a, 
2003b). 

Recently, Magne et al. (2005, hereinafter referred to as MAHR) showed that AH's 
theory gives the same damping of incident waves as the Green function solution of Pihl, 
Mei & Hancock (2002), applied to any two dimensional topography, random or not. 
Investigating the applicability limits of the scattering term of AH, MAHR also performed 
numerical calculations, comparing AH's theory to the accurate matched-boundary model 
of Rey (1992) that uses a decomposition of the bottom in a series of steps, including 
evanescent modes. The numerical results show that AH's theory is generally limited by 
the relative bottom amplitude 77 = max.{h)/H rather than the bottom slope. In particular, 
AH's theory predicts accurate reflections, with a relative error of order 77, even for isolated 
steps that have an infinite slope (MAHR). 

The resulting expression of the Bragg scattering term is consistent with results for scat- 
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tering of acoustic and electromagnetic waves obtained by the small perturbation method, 
valid in the limit of small k max(/i) where k is the wavenumber of the propagating waves 
(Raylcigh 1896, see Elfouhaily & Guerin 2004 for a review of this and other approxima- 
tions). Since there is no scattering for kH » 1, as the waves do not 'feel' the bottom, 
the small parameter n = max(/i)/if may be used in the context of surface gravity waves, 
instead of the more general fcmax(/i). For n <C 1, the scattering strength is thus entirely 
determined by the bottom elevation variance spectrum at the bottom scales resonant 
with the incident waves. 

Based on these results, Mei's (1985) theory should yield the same reflection coefficient 
as AH's theory in the limit of small bottom amplitudes. Yet, AH predict that the wave 
amplitude in 2D would decay monotonically, which is not compatible with the oscillatory 
nature of Mei's theory for large detunings from resonance. Further, outside of the surf 
zone and the associated multiple bar systems, the application of AH's theory is most 
relevant in areas where the bottom topography changes significantly on the scale of the 
wavelengths of swells. This often corresponds, over sand, to the presence of sandwaves. 
These sandwaves are generated by currents, and particularly by tidal currents (e.g. Dal- 
rymple Knight & Lambiase 1978; Idier, Erhold & Garlan 2002). It is thus logical to 
include the effects of currents in any theory for wave scattering over a random bottom. 
Kirby (1988) developed such a theory for monochromatic waves over a sinusoidal bottom 
and a slowly varying mean current, extending Mei's (1985) work. The geometry of the 
resonant wavenumbers is modified in that case, with with incident and reflected waves 
having the same absolute frequency, but different wavenumber magnitudes if incident 
and reflected waves propagate at different angles relative to the current direction. Kirby 
(1988) also considered the short scale fluctuations of the current, due to the sinusoidal 
bottom, that may be interpreted as a separate scattering mechanism, and generalized fur- 
ther to any irrotational current fluctuations, leading to results similar to those obtained 
for gravity-capillary waves by Bal & Chou (2002). 

The present paper thus deals with these two questions. An extension of AH's theory 
for surface gravity wave scattering in the presence of irrotational currents with uniform 
mean velocities is provided in § 2, and the differences between this theory and those 
of Mei (1985) and Kirby (1988) are discussed in detail in § 3. Expected oceanographic 
effects of scattering in the presence of a current are investigated in § 4, using a spectral 
phase-averaged numerical model, predicting the evolution of the wave action spectrum, 
and detailed measurements of the topography in the southern North Sea. Conclusions 
follow in § 5. 

2. Theory 

2.1. General formulation 

The variation in the action spectral density due to wave-bottom scattering is derived 
following the method of AH, now including the effect of a uniform mean current. The 
method is identical to that of Kirby (1988) with the difference that an equation for the 
spectral wave action is sought instead of one for the wave amplitudes. Thus intermediate 
results are identical to those of Kirby (1988). Since the wave action is a quadratic function 
in the wave amplitude, we will naturally consider the wave potential up to second order 
in the normalized bottom amplitude rj, in order to have all wave action terms to order 
rj 2 . The only important terms in this type of calculation are the 'secular terms', i.e. 
the harmonic oscillator solutions for the wave potential forced at resonance, with an 
amplitude that grows unbounded in time. We shall thus obtain a rate of change of the 
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Figure 1. Definition sketch of the mean water depth H, and reiative bottom eievation h, for 
one particular case of a sinusoidal bottom investigated in § 3. 



action from the equality of all the secular terms. The particularity of the random wave 
approach is also that we will consider all possible couplings between wave components, 
and not just two wave trains. With random waves, secularity is limited to a sub-space of 
the wavenumber plane that generally has a zero measure. Thus the near-resonant terms, 
once integrated across the resonant singularity, are the ones that provide the secular terms 
for random waves. This integration assumes that the spectral properties are continuous, 
a real theoretical problem for nonlinear wave- wave interactions (e.g. Benney & Saffmann 
1966, Onorato et al. 2004). Here we shall see that the only relevant condition is that the 
bottom spectrum be continuous, at least in one dimension. This is obviously satisfied 
by any real topography, since a truly infinite sinusoidal bottom of wavelength L, with 
an infinite spectral density at the wavenumber 2n/L, is not to be found, even in the 
laboratory. 

We consider weakly nonlinear random waves propagating over an irregular bottom 
with a constant mean depth H and mean current U, and random topography h(x), with 
x the horizontal position vector, so that the bottom elevation is given by z = — H + h(x) 
where z is the elevation relative to the still water level. The bottom undulations cause 
a stationary random small-scale current fluctuation (u(x, z), w(x, z)) deriving from a 
potential <p c . The free surface is at z — <^(x, t) — C(x, t). Extension to mean current and 
mean depth variations on a large scale follows from a standard two-scale approximation, 
identical to that of by Kirby (1988). This is not included in the present derivation for 
the sake of simplicity. 

The maximum surface slope is characterized by e and we shall assume that e 3 <C i] 2 
so that the bottom scattering contributions to the wave action to order rf are much 
larger than the resonant non-linear four wave interactions (Hasselmann 1962) that shall 
be neglected. Such interactions could also be allowed in the present calculation providing 
an additional source of scattering with the known form due to cubic non-linearities. For 
shallow water waves {kH << 1) a stricter inequality is needed to prevent triad wave- 
wave interactions to enter the action evolution equation at the same order as bottom 
scattering. 

The solution is obtained in a frame of reference moving with the mean current vector 
U , which has the advantage of removing the convective terms due to the mean current 
velocity. The corresponding transformation of the horizontal coordinates is x 1 = x + Ut, 
where x and x' are the coordinates in the moving and fixed frames, respectively. As a 
result of this transformation, the bottom is moving, and the bottom boundary condition 
for the velocity potential is modified. The governing equations consist of Laplace's equa- 
tion for the velocity potential, which includes both wave and current motions, the bottom 
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kinematic boundary conditions, and Bernoulli's equation with the free surface kinematic 
boundary condition. Assuming that the atmospheric pressure is zero for simplicity, and 
neglecting surface tension, one has 

V 2 (f>+—^=0 for -H + h^z^C, (2-1) 



dz 2 



dz dt 



dh 

V0-V/i at z = -H + h, (2.2) 
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f +*C=4lW| 2 -i(^ + c(t) at , = C . (2.3) 

| = | + V0.VC at z = C, (2.4) 

with c(t) a function of time only, to be determined. The symbol V represents the usual 
gradient operator restricted to the two horizontal dimensions. The latter two equations 
may be combined to remove the linear part in £. Taking <9 lj2.3|l /(9t +g(J23J, yields, 



d 2 (t> dej) _ , 8C, d 2 cj) ( dC, d 

w+9 -=gV<t>.Vt--—-ll+-- 



dV<p d<f) d 2 <\> 



c'(t), 



dt dz dtdz 

at z = C (2.5) 



Following Hasselmann (1962), we approximate h and <f> with discrete sums over Fourier 
components, and take the limit to continuous integrals after deriving expressions for the 
evolution of the phase-averaged wave action. We look for a velocity potential solution in 
the form 

0(z, M ) = £ 8i(z, Tt)^**-^ = £ g&W ^^f^ *-' + ■ ■ • - ( 2 -6) 

k.s k,s 

where a is the radian frequencies in the moving frame, k is the surface wavenumber, 
with magnitude fc, and s is a sign index equal to 1 or — 1. In the moving frame of 
reference, s = 1 for wave components that propagate in the direction of the vector fc, 
and s = — 1 for components that propagate in the opposite direction. Thus the radian 
frequency in the fixed frame is u> — a + sk • U. The amplitudes are slowly modulated 
in time, with a slowness defined by the small parameter 7. Because is a real quantity 

we also have = $1^., where the overbar denotes the complex conjugate. Thus the 
double decomposition made in Ij2.6|l in wavenumber k and propagation direction + or 
— replaces a more general decomposition in wavenumber and frequency that would be 
necessary if nonlinear dispersive effects were included. Here the frequency a is always 
related to k via the linear dispersion relation. 

In the alternative decomposition with amplitudes ^ that contain the fast time vari- 
ation, only the part of the solution that has the vertical structure of Airy waves has 
been given explicitly. The other part, represented by '. . .', will be found to be negligible 
for small bottom amplitudes. Our rather archaic use of the s index to distinguish the 
wave propagation direction is preferred here to the more modern use of the Hamilto- 
nian variables that combine elevation and potential at the free surface, widely used for 
wave- wave interaction studies (e.g. Janssen 2004). The complexity of the Hamiltonian 
variables appears unnecessary for the linear waves considered here. 

Expanding the bottom boundary condition and wave potential in powers of 77 = 
max(h)/H, 

= 0o + 0i +^2 + ..., (2.7) 
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where each term fa is of order rf. The boundary conditions l|2.5|) and Ij2.2(l are expressed 
at z = and z = —H, respectively, using Taylor series of <f> about z = —H and z = 0. 

Unless stated otherwise, these potential amplitudes will be random variables. Since 
we are solving for </> seeking an equation the for wave action N, we must relate N to 
4>. Accurate to second order in e and r\ (see Andrews & Mclntyre 1978 for the general 
expression of N) we have N = E/a for a monochromatic wave of surface elevation 
variance E and intrinsic frequency a, in which, following the common usage in non- 
accelerated reference frames, the gravity g is left out, so that the action has units of 
meters squared times second. For general waves, the variance E may be written as 

E(t) = ((Co + Ci + C 2 + • ■ -) 2 ) = (Co 2 + 2CoCi + (Ci + 2C0C2) +..•), (2.8) 

where (•) denotes an average over flow realizations, and Q is the surface elevation solution 
of order ?y l , and terms of like order in 77 have been grouped. Each of these terms may be 
expanded in this form 

(Co 2 ) = E fef = 2$>oVo:- fc (2.9) 

fc,s k 

For free wave components, the elevation amplitude is proportional to the velocity poten- 
tial amplitude 

Z',H = isv*U/9 (2-10) 
so that the elevation co- variances are proportional to the co-variances F^ k of the surface 
velocity potential, 

F*U = "W";. A: + (2-H) 

The contribution of the complex conjugate pairs of components (fe,+) and (— k, — ) 
are combined in l|2.11|l so that the covariance i^* k correspond to that of all waves 
with wavenumber magnitude k propagating in the direction of k. In the limit of small 
wavenumber separation, a continuous slowly-varying cross-spectrum can be defined (e.g. 
Priestley 1981, ch.ll; see also AH), 

F*Jk)= lim A fV, ■ (2.12) 
' |Afc|->o Ak x Ak y 

The definition of all spectral densities are chosen so that the integral over the entire 
wavenumber plane yields the total covariance of fa and <fij . 

Finally, Nij(k) is defined as the (i + j) th order depth-integrated wave action contribu- 
tion from correlation between i th and j th order components with wavenumber k. From 



lf2~H|) and lf2~TU|l one has, 



N itj {k) = —F^(k)t£inh(kH). (2.13) 



The spectral wave action is thus, 

oc 00 i 

N(k) = J2^i(k) =^^% f (fe). (2.14) 

i=a i=0 j=0 

Defining Gi as the amplitude of the Fourier component of wavenumber I, the bottom 
elevation is given by 

h(x) =E G ' e "' [x+C/tI ' ( 2 - 15 ) 
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with a summation on the entire wavenumber plane. Because h is real, G_j = Gj. The 
bottom elevation spectrum in discrete form is given by = (GiG-i) and in continuous 
form by 

F B (l)= lim -A- , (2.16) 

^ ' |A/|->0 Al X Aly 

and verifies, 

/oc />oc i rL/2 rL/2 

/ F B (I)dI x dZ» = lim / / h 2 (a:, y)dzdy (2.17) 
-ooJ-oo L^oo L J^ L / 2 J-L/2 

Now that the scene is set, we shall solve for the velocity potential <f> in the frame of 
reference moving with the mean current, and use (|2.13(l to estimate the action spectra 
density at each successive order. In the course of this calculation, <f> will appear as the sum 
of many terms, some of which are secular (these are the 'resonant terms' in Hasselmann's 
terminology), i.e. with growing amplitudes in time. Most importantly among these are 
those that lead to resonant terms in N. All other terms are bounded in time and thus 
do not contribute to the long-term evolution of the wave spectrum, i.e. on the scale of 
several wave periods, and shall be neglected (see Hasselmann 1962). 

2.2. Zeroth-order solution 

In the moving frame of reference, the governing equations for <po are identical to those 
in the fixed frame in the absence of current. The solution is thus 

E cosh(ib(z + H)) ^. a _ mt] 
t-^ cosh(kH) 

k,s 

where the intrinsic frequency a is the positive root of the linear dispersion relation, 

a 2 = gktanh(kH). (2.19) 

2.3. First-order solution 

Surface non-linearity becomes relevant at first order due to a coupling between the ze- 
roth order solution and current-induced first order terms. Including all powers of rj, the 
expansion of the surface boundary condition to order e 2 gives, at z = 0, 

d</> | d</> _ #V ^d 2 <f> d( d 2 (/> | v ^ / dV<A dcj) d 2 <j> 

dt 2 dz dt 2 dz dz 2 dt dzdt \ dt J dz dtdz 

+c[(t) + 0(e 3 ) (2.20) 
The equations at order r\ are 

V 2 0i + — % = for - H ^ z sC 0, (2.21) 

dz z 

^ = -h*& + V*o-Vh+% at z = -H, (2.22) 
dz dz 1 at 

and, at z = 0, expansion of l|2.20l) to first order in r/ yields, 

i ii 



in 



^ +9 ^z~= 9 ( v v ^- v C^Ci^J+. 9 ( v V0 1 .vCo-Co^J — 

IV V VI VII VIII 

_V0 O • 9V(t>1 (— + —) ^_ 2 ^l^i_ Cl ^L _ Cf) f^L +NLl 
dt \ dz dt J dtdz dz dtdz dt 2 dz dt 2 dz 
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(2.23) 

where the terms NL\, not written explicitly (see Hasselmann 1962 eq. 1.11-1.12), are 
quadratic products of the zeroth-order solution. Since no gravity waves satisfy both 
a = G\ ± (i2 and k = k\ ± fc2, NL\ forces a non-resonant wave solution t^" 1 that will be 
neglected because it does not modify our second order wave action balance, thanks to the 
choice e < rj. The spatially uniform term c' 2 (i) has been incorporated into NL\ and is also 
of second order in the wave slope, and does not lead to resonances. That term, omitted 
by Hasselmann (1962), is responsible for generating microseisms (e.g. Longuet-Higgins 
1950). 

The hrst-order system of equations is non-linear due to the surface boundary condition 
H2.23[l . However, all the right hand side terms of (|2.23(1 are of order er/tpo, and thus 
negligible, provided that 4>i is of order r](j)Q. Without dh/dt in H2.22(l this would be the 
case, since the other forcing terms are all proportional to r]4>o- However, as suggested 
by anonymous reviewers, dh/dt introduces an external forcing. We thus first give the 
solution (</>i c , Cic) forced by dh/dt only, in the right hand side of l|2.22|l . This solution is 
physically identical to mean current perturbation caused by the bottom topography and 
given by Kirby (1988, his eq. 2.9) for a sinusoidal bottom. With a more general bottom, 
it is 

<j)i c = i Vf7 • I— - {# cosh [l(z + H)] + ai sinh [l(z + H)]}^ x+Ut \ (2.24) 
^ lai 



where 



and 



(U • I)' 2 

ai = i - t&nh(lh), (2.25) 

A = 1 - tanh(l/t) ( U 'P . (2.26) 

gi 

The corresponding surface elevation oscillations, given by (|2.3p . are second order in the 
Froude number Fr — U / '(gh) 1 / 2 , and 180° out of phase with the bottom oscillations for 
slow currents when a < (Kirby 1988, eq. 2.10), 

Ci c = E {U ' l) lnL eiHx+m) - ^ 

' gai cosh(m) 

From l|2.24|l . the following expression are derived, 

ihc{z = o) = iJ2 u • h — %m ciHx+ut) > ( 2 - 28 ) 

*— J lai cosh(m) 

?pl {z = o) = % = iV (U ■ if (2.29) 

dz v ' dt ' gaicosh(lh) v ' 

These shall be particularly useful for plugging into l|2.23ll . 

We can now obtain the general solution to our equations (|2.21|) - (|2.23|l by the follow- 
ing superposition of the previous solution with free and bound (i.e. non-resonant) wave 
components, with amplitudes ^\ k and respectively, 



E 



cosh [k(z + H)} sinhj^z + Hl 8ijS . 



cosh(kH) L '" y ' cosh(kH) 



c ik ' x , (2.30) 
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where the last two terms corresponds to the solution to the forcing by all the right hand 
side terms except for dh/dt. Because <j>i c and Cic are the only terms that may be larger 
than eijipo, all others are neglected in the right-hand side of Q2.23p . 
Substitution of (|2.3U|) in the bottom boundary condition 12.22fl yields 



cosh(fctf) 1,fcW ^-fcosh(k'H) 

Replacing now l|2.3U[) in the surface boundary condition H2.23[l . yields an equation for 
$1 k . Using id = a + k • U and u>' = a' + k' • U, it writes 



alt 2 



A = E AfS ( fc ' k')^ Q yG k _ k A h - u -^ ']*, (2.32) 

/ 1 1 



with 



M s (k,k') = \gk-[k-U- su;'} 2 tanh(kH)} - ' k + M s cl {k, k') (2.33) 

where is given by all the right-hand side terms in (|2.23() and thus corresponds 
to the scattering induced by current and current-induced surface elevation variations. 
Anticipating resonance, we only give the form of = for a = a' — si • U , with 

(d) 

(a) (b) „ (c) 



_ {sg 2 U ■ l(a'l .k + al.k')-(U. if [g 2 k ■ k> - oo\oo' + (U ■ I) 2 )}} 

Mc(fe ' fej " lg^aicosh(lh) 

(2.34) 

in which the term (a) is given by the term (II) in H2.20JI . (b) is given by (III) and (IV), (c) 
is given by (I), and (d) is given by (V)-(VIII). Because we are first solving the problem 
to order 77, it is natural that our solution is a linear superposition of the solutions found 
by Kirby (1988) for a single bottom component. Indeed, M c (k, k') — — 4ui£l c /D, with il c 
the interaction coefficient of Kirby (1988, eq. 4.22b) and D his bottom amplitude, here 
Gi = LD/2. 

The solution to the forced harmonic oscillator equation (|2.32|1 is 

*;,*(*) =E MS ( fe ' fe >o, fe ' G fe-fe'/i(< M ■ u ~ sa ' ;t) > (2 - 35) 
k' 

where I = k — k' , and the function f\ is defined in Appendix A. 

2.3.1. First order action 

The lowest order perturbation of the wave action by scattering involves the order n 
covariances 

F? <0 , h + F * lifc = 4Re ((^Z-k)) , (2-36) 
with Re denoting the real part. Including only the secular terms, we get 



Y,M + (k,k?)($+^_ h G h _M(a,l.U-</;t)<& 
L k' 



(2.37) 



Although this term was assumed to be zero in AH, it is not zero for sinusoidal bottoms 
with partially standing waves, and may become significant at resonance due to the func- 
tion /1. In uniform conditions, the time evolution of the wave field requires that the 
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non-stationarity must come into play. Thus 7 « 77 and the non-stationary term is given 
by AH (their appendix D), 

d[N%>(k) + NS\(k)] dN (k) 

dt — ^ — - w^- (2 ' 38) 

In order to simplify the discussion, we shall briefly assume that there is no current and 
that the waves are unidirectional. In that case, fc' = — fc and M(k, fc') = —gk 2 / cosh 2 (kH). 
Replacing H2.37[l in (|2.13(1 and combining it with (|2.38|1 yields the action balance 



dN , k = d_ 
dt dt 



A tanh(A:i/ ) (F* 0>k + F* hk ) 



4fc 2 cr 

(2.39) 

with Im denoting the imaginary part. 

For directionally spread random waves, with a current, and a real bottom (e.g. random 
or consisting of a finite series of sinusoidal bars), the evaluation of (|2.37() is not simple. 
First of all, resonant terms given by /1 only occur for a' = a + si • U , that is uj = lu' . 
Using N(k) — N (k) [1 + 0(r])] and taking the limit to continuous surface and bottom 
spectra yields 

with the mixed surface bottom bispectrum Z defined by 

$ i"fe*7 h 'G-k-k> 

with k = k(cos9,sm9) and fc' = k(cos9' ,sm6'). Z is similar to a classical bispectrum 
(e.g. Herbers et al. 2003) with one surface wave amplitude replaced by a bottom ampli- 
tude, and a similar expression is found for a non-zero current. The action balance l|2.40|l 
is generally not closed, and requires a knowledge of the wave phases that are not available 
in a phase-averaged model. The same type of coupling, although due to the large scale 
topography, also occurs in the stochastic equations for non-linear wave evolution derived 
by Janssen, Herbers & Battjes (2006). 

The contribution of the mixed bispectrum will thus be evaluated below, in order to 
investigate in which cases it may be neglected or parameterized. It is expected that S\ is 
generally negligible because MAHR have neglected Si , and still found a good agreement 
of the second order action balance with exact numerical solutions for the wave amplitude 
reflection coefficient. 

2.3.2. Second order action 

From the expansion l)2.14[l . the second order action is A^fc) = Nx t i(k) + No^(k) + 
A^2,o(fc)- The first term can be estimated from cj>x, using the covariance of the velocity 
potential amplitudes l|2.11(l . 

F?A,k = ZiK^i-k)- ( 2 - 42 ) 
Using H2.35fl . (|2.42|l can be re- written as 

2 



F 



O.fc' 



' (l^fe-fe'G-fe+fc' 



2, 



^ - 2"^ |M+(fc,fc')| 2 ' ' ^''W'^ 1 ; \f 1 (a,l-U-a f ;t)\ A Ak', 

(2.43) 



Afc ^V 1 ' Afc' Afc 

k' 
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Taking the limit of (|2~I3|) when Afc -> 0, 

/OO />OC 
/ \M+{k, k')\ 2 F^ 1 {k')F B {k - k') |/i(cr, I ■ U - a'; t)f dk' x dk' y . 
-OO J — OO 

(2.44) 

Due to the singularity in fi , and assuming that the rest of the integrand can be approx- 
imated by an anlytical function in the neighbourhood of the singularity u>' = to, which 
requires both bottom and surface elevation spectra to be continuous, the integral can be 
evaluated by using 

nt 

I • U - a'; t)fi(a, 1-U + a'; *)) = — [8 (o' - (a + I ■ U)) + 0(1)] . (2.45) 

S is the one-dimension Dirac distribution, infinite where the argument is zero, and such 
that J 5(x)A(x)dx = A(0) for any continuous function A. In order to remove that singu- 
larity, the argument of S maye be re- written as ui'—u, making explicit all the dependencies 
on k' . Evaluation of the 5 function is then performed by changing integration variables 
(k' x ,k' y ) are changed to (to', 9'), with a Jacobian k'dk'/dw' = k' 2 / {k'C' g + k' -U). We thus 
have 

ret f 2n f o k'F B (k — k'\ 

F? 1 (t,k) = ^ j( ^ |M+(fc, fc')| 2 F*(fc') c , f$ . ^ V ' - u>)du/M + 0(1). 

(2.46) 

When uj = u)\ the integrand simplifies. M s (k, k') is equal to M(fc, fe'), defined by 

ok • k' 

M(k, k') = y — — + M c (k, k') = M b (k, k') + M c {k, k% (2.47) 

cosh(fcii) cosh(/c H) 

with M c = M+ given by i|2.34[l . Using the (|2.13l) relation between velocity potential and 
action, and evaluating the integral over u> , one obtains 



2 



Again we note the correspondance with the theory of Kirby (1988, eq. 4.21). Specifically, 
one has M(k, k') = —4ui£l c / D, with il c being Kirby's interaction coefficient. 

2.4. Second order potential and corresponding terms in N2 

In order to estimate the other two terms that contribute to N2, the second order potential 
4>2 must be obtained. It is a solution of 

d 2 d> 

V 2 2 + -rrv =0 for -H ^z^O, (2.49) 
oz z 

d(f> 2 ,d 2 <px h 2 d 3 (j> ryru d< t>0\ . xi hem 

57 = " /l 9^-Y^ + V01 ' v/l + v(/l 97 ) - v/l at Z = - H > (2 - 50) 

that simplifies because odd vertical derivatives of 4>o are zero at z = —H, 

^l = _ h ^l + V<j)x . Vh at Z = _ H , (2.51) 
oz oz 1 

and 

d 2 cf> 2 dfo d$ 



dt 2 9 dz * dt 



= i 2sa—^-e i ^- !D - su ^ + I - VIII + NL 2 at z = 0. (2.52) 



12 R. Magne and F. Ardhuin 

The terms I- VIII are identical to those in Ij2.23|) with tfio, Co, 4>i, Ci replaced by (<fii — 4>x c ), 
(£1 — Cic), 4*2 and C2, respectively. All other non-linear terms have been grouped in NL,2- 
In order to yield contributions to the second order action A^o, terms must correlate with 
4>o to give second-order terms in r\ with non-zero means. For zeroth order components 
with random phases, inspection shows that NL2 do not contribute to A^.o and will thus 
be neglected. 

The solution <p2 is given by the following form, 



E 

fe.s 



cosh(fcp + H )) ^ s siiih(fc(z + g)) S| ' 

cosh(kH) ^ k{t) + co S h(kH) 2 > k[t \ 



c lk ' x . (2.53) 



The non-stationarity term </>!} s leads to the action evolution term (|2.38(1 . now assuming 
7 ss rj 2 . Following the method used at first order, substitution of (|2.53(l in the bottom 
boundary condition 12.51J1 leads to, 

•SCO = -E i T i ^T^.,WG,-.e»- . (2.54) 

After calculations detailed in Appendix B, <f>2 yields the following contribution to the 
wave action, 

N 2 , (k) + N , 2 (k) = ~ I M 2 (fc, k')F B (k k' fy — *- -, v d9' + 0(1), 

(2.55) 

in which a' = a - I ■ U, a' 2 = gk'tanh(kH), and C' g = t/(l/2 + k'H/ smk(2k'H))/k'. 

2.5. Action and momentum balances 

We shall neglect the first order action contribution iVi given by H2.40JI . The solvability 
condition imposed on the action spectrum is that A^ remains an order rj 2 smaller than 
iVo for all times. Thus all secular terms of order ?/ 2 must cancel. Combining 
and H2.55fl gives 

d7V (fc) , n f 2n nr2n „ ^„ B ,u ,,, N o(k')-N Q (k) 



Since N2 and Ni remain small, N(k) — Ao(fc) [l + 0(?? 2 )] , and one has, 

dTV(fc) 



(2.56) 



dt 

with the spectral action source term, 

,27t k' 2 M 2 (k,k') 



'bscat 



(fc), (2.57) 



Si 



7T f zn V l M z (h k'\ 
bsca t (fc) = £ / 7(un. :l TA ^^ - - N W d0 '> ( 2 - 58 ) 

2 Jo o-a' (k'C'g + k'-U) 

where a' = cr+l-U and k = k'+l. This interaction rule was already given by Kirby (1988). 
The only waves that can interact share the same absolute frequency ui = a + k ■ U = 
a' + k' ■ U. For a given k and without current, the resonant k' and I lie on circles in the 
wavenumber plane (see AH). The current slightly modifies this geometric property. For 
U << C g the circles become ellipses (Appendix C). 

For a given value of u>, one may obtain the source term integrated over all directions, 

f 2n dk 

SWatM = / kS hscat {k)—d9 (2.59) 
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fte f2n n k 2 k' 2 M 2 (k,k')F B (k-k') 



J 2 aa' (k'C' + k' -U) [kC g + k-U) 



[N(k') - N(k)} d9'd6 



2n r 2n i\j2fu U'\ 



7t M 2 (k,k')F B (k - k') 



k 2 N(uj,9') k' 2 N{u,9) 



kC g + k-U k'C' g + k'-U 



de'cie. 



This expression is anti-symmetric, multiplied by -1 when 9 and 9' are exchanged. Thus 
"Sbscatfw) is a substraction of two equal terms, so that for any bottom and wave spectra 
5bscat(^) = 0. In other words, the 'source term' is rather an 'exchange term', and con- 
serves the wave action at each absolute frequency. This conservation is consistent with 
the general wave action conservation theorem proved by Andrews & Mclntyre (1978), 
which states that there is no flux of action through an unperturbed boundary (here the 
bottom). It also appears that ui and 9 are natural spectral coordinates in which the 
scattering source term takes a symmetric form. Finally, we may consider the equilibrium 
spectra that satisfy Sbsc&t{k) — for all k. Without current, an equilibrium exists when 
either N(ui, 9) or N(k) is isotropic. With current, the scattering term is uniformly zero 
if and only if the spectral densities in fc-space, N(k), are uniform along the curves of 
constant to. 

The source term SWat may also be re-written in a form corresponding to that in AH, 
which now appears much less elegant, 

S bs cat(fc)= / K(k,k , ,H)F B (k - k')[N(k') - N(k)}d9', (2.60) 
Jo 

with 

„ ,, , nk' 2 M 2 (k,k') Anakk' 3 cos 2 (9 - 0')[1 + O(Fr)] 

K(k, k , H ) — 



2aa' (k'C' g + k' • U) sinh(2fc J ff) [2k' H + sinh(2fc'ff) (1 + 2k' ■ Ufa')] ' 

(2.61) 

One may wonder how large is the current-induced scattering represented by M c , our 
eq. Q2.34p . compared to the bottom- induced scattering represented by Mb. Since a' = 
er + (U • I), the (a) and (b) terms in the numerator M c almost cancel for small Froude 
numbers, and the (a)+(b) part is of order Fr 2 . Thus M c is generally an order Fr 2 
smaller than Mb. For k and k' in opposite directions (i.e. back-scattering), the (a)+(b) 
part is even smaller, of order g 2 (U • l) 3 l ■ k, and exactly zero in the long wave limit 
IH <C 1. Thus, for back-scattering, the numerator in M c is itself of the order of (c), i.e. 
g 2 (U • l) 2 k ■ k' . Interestingly (c) formally comes from the modulations of the surface 
elevation £ic so that the 0(Fr 2 ) elevation modulation is at least as important as the 
O(Fr) current modulation for this back-scattering situation. In that case, M c is of the 
order of Mb cosh(fcff) cosh(k' H)(U • l) 2 /[glai cosh(lH)]. The relative magnitudes of Mb 
and M c thus depend on Fr(l) = {U -l) /{glt&nh(lH)\ 1 / 2 that appears in (U • I) 2 / \gla{) = 
Fr 2 (l)/[Fr 2 (l) - 1]. This /-scale Froude number may be formally close to 1, and thus 
M c may be larger than Mb. However, scattering is limited by blocking as no scattered 
waves can propagate when Cg' < U ■ k'/k'. In the long wave limit, Fr(l) — Fr and 
for (1 — Fr) <C 1, one has M c > Mb- For oblique scattering, the (a)+(b) term may 
dominate the numerator of M c and the situation is more complex. Nevertheless, for 
Froude numbers typical of continental shelf situations, say < Fr < 0.4, M c may be 
neglected in most situations since its 0(Fr 2 ) correction corresponds to only a few percent 
of the reflection. Obvious exceptions are cases in which Mb is zero, such as when k and 
k' are perpendicular. 

Finally, we may also write the evolution equation for the wave pseudo-momentum 
= p w g J kN(k)dk (see Andrews & Mclntyre 1978), where p w is the density of 
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sea water. Introducing now the slow medium and wave field variations given by Kirby 
(1988), that do not interfere with the scattering process, except by probably reducing 
the surface-bottom bispectrum Z , one obtains an extension of the equation of Phillips 
(1977) 



~dT dx^ [{Up + Gsf))Ma ] ~~ Ta ~ M ? dx^ ~ 1^^2kDdx~ a 



with the dummy indices a and (3 denoting dummy horizontal components, and the scat- 
tering stress vector, 



r bscat = - Pw g J fcS bscat dfc. (2.63) 

This stress has dimensions of force per unit area, and corresponds to a force equal to 
the the divergence of the wave pseudo-momentum flux. Based on the results of Longuet- 
Higgins (1967) and Hara & Mei (1987), this force does not contribute to the mean flow 
equilibrium with a balance of the radiation stresses divergence by long waves (or wave set- 
up in stationary conditions), contrary to the initial proposition of Mei (1985). This force 
is thus a net flux of momentum through the bottom, arising from a correlation between 
the non-hydrostatic bottom pressure and the bottom slope. That force is likely related to 
the pressure under partial standing waves locked in phase with the bottom undulations. 
Although the part M c of the coupling coefficient M given by l|2.47|l is formally due to 
scattering by the current modulations V0i c , and associated surface fluctuations £i c , it 
should be noted that these motions and related pressures are correlated with the bottom 
slope in the same way as the part represented by M^. Thus both terms contribute to this 
force r bscat which acts on the bottom and not on the mean flow. 



3. Wave scattering in two dimensions 

Before considering the full complexity of the 3D wave-bottom scattering in the presence 
of a current, we first examine the behaviour of the source term in the case of 2D sinusoidal 
seabeds. Although the bottom spectrum is not continuous along the y-axis, continuity 
in x is sufficient for the use of (|2.45l) and the source term can be applied, after proper 
transformation to remove these singularities. MAHR have investigated the applicability 
limits of the source term with [7 = 0. They proved that for small bottom amplitudes the 
source term yields accurate reflection estimates, even for localized scatterers, and verified 
this with test cases. It is thus expected that this also holds for U ^ 0. 

3.1. Wave evolution equation in 2D 

We consider here a steady wave field in two dimension with incident and reflected waves 
propagating along the a;-axis. We shall consider in particular the case of m sinusoidal 
bars of amplitude b and height 2b, with a wavelength 2tz/Iq. The bottom elevation is thus 

h(x) = bsm(mlox) for < x < L (3-1) 
h(x) = otherwise. 

Such a bottom is shown in figure 1 for m = 4. This form is identical to that of the bottom 
profile chosen by Kirby (1988) but differs, for < x < L, by a n/2 phase shift from the 
bottom profile chosen by Mei (1985). The bottom spectrum is of the form 



F B (l x ,l v ) = F B2D (l x )S(l y ), 



(3.2) 
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and for the particular bottom given by 1)3. l|l . 



F 



B2D 



1 

2n 



h(x)e~ llx dx 



2b 2 l 2 sin 2 (lL/2) 



7tL 



(11 



l 2 f 



(3.3) 



F 



B2D 



(3.4) 



with 

mb 2 _ b 2 L 
"44" ~ ~8n' 

Note that this is a double-sided spectrum, with only half of the bottom variance contained 
in the range l x > 0. For a generic bottom, for which h(x) does not go to zero at infinity, 
the spectrum is obtained using standard spectral analysis methods, for example, from 
the Fourier transform of the bottom auto-covariance function (see MAHR) . In that case 
pB2D j g e q U i va i en t to a Wigner distribution (see e.g. Ryzhik et al. 1996). 

First, replacing 1)3. 2|) in l|2.57[) removes the angular integral in the source term. Taking 
k = (k x , k v ), we have l y = k y 



k' y = ksm9 — k'sm.6', thus dl y — —k'cosO'dO', and 



5*bscat (fe; X) 



nk'APjk^^F^ik-k') 
2crcr'|cos6»'| (k'C'g + k' • U) 



[N{k') - N{k)\ . 



(3.5) 



Second, assuming now that waves propagate only along the cc-axis, the wave spectral 
densities are of the form 



N(k x ,k y ) = N(k x ,ky)S(k y ) = N 2U (k)S(0 - 9 )/k 



(3.6) 



with da = for k x > and 6q = n for k x < 0. Integrating over 8 removes the singularities 
on k y , and assuming a steady state one obtains 



k x n 



u x 



dN 



2D 



dx 



(k x ,x) — Sbscat i^X) X) , 



with 



n2D 



bscat 



(k x , x) 



nk'M 2 (k,k')F B2U (k x -k' x ) 
2aa' (k' x C' g + k' x U x ) 



N 2D (k' x ,x)-N 2D (k x ,x)} 



(3.7) 



(3.8) 



Although the present theory is formulated for random waves, there is no possible 
coupling between waves of different frequencies. Mathematically, it is possible to take 



the limit to an infinitely narrow wave spectrum, such that, N 2D 
N'(x)8{lu' with k 0x > and k' 0x < 0. Using duj/dk = C g 4 

evolution equation is, omitting the subscripts on k and k' , 



k x n 

T° 9 



u x 



dN 
dx 

nM 2 (k,k')F B2D (k x 
W 



kN' 



(k,x) = N(x)S{uj-uj a ) + 
-k x U x / \k x \, the resulting 



k'N 



kC q + k x U x 



k'C'g + k x U x 



(3.9) 

with a similar equation for N' obtained by exchanging C g and C' g , and k' and k, from 
which it is easy to verify that the total action is conserved. 

The stationary evolution equation (|3.7|l only couples two wave components N(k) and 
N(k'). For a uniform mean depth H, and uniform bottom spectrum F B , as considered 
here, we thus have a linear system of two differential equations, that may be written in 
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matrix form for any k > 0, 
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N(k) 



qQ 



N(k) 
N(k') 



with 



dx \ N(k') 

nM 2 {k,k')F B2D (l) 



(3.10) 



(3.11) 



2aa'CgCg' 

Defining I = k x — k' x , the action advection velocities V = C' g +k' x U x and V = C g +k x ll : 
the terms of the non-dimensional matrix Q are given by 



(0)2,1 



c 9 c> 

V 2 
V' 2 



and (Q)i,2 
and (0)2,2 



vv ' 

W ' 



(3.12) 



where (0)j,j is the r row and j column term of Q. The general solution is thus 



N(k,x) 
N(k',x) 



Qx ( N(k,0) 
N(k',0) 



(3.13) 



The matrix exponential is classically the infinite series X^^Lo (?0) in which matrix 
multiplications are used. The reflection coefficient for the wave action is found using the 
boundary condition expressing the absence of incoming waves from beyond the bars, 
N(k',L) = 0, giving, 



N(k',0) 
N(k,0) 



(3.14) 



A reflection coefficient for the modulus of the wave amplitude predicted by the source 
term is thus, 



Rs = 



a'N(-k',0) 



aN(k,0) 



1/2 



{°'(« Q %i/[<'(' 9QL )w]Y' a - ( 3 - 15 ) 



The spatial variation of the amplitudes may be linear, oscillatory, or exponential, de- 
pending on whether the determinant of Q, is zero, negative or positive, respectively. 
That determinant is C 2 C' 2 (V - V)(V' 2 + 3VV + W^/V^V' 3 , which is always of the 
sign ofV' — V. 

3.2. Analytical solution for U = 
In the absence of a mean current, k! = —k, and 



(<?)i,i = (Q)i 



-(0)2,1 = (0)i 



(3.16) 



Thus Q 2 = so that its exponential is only the sum of two terms, e qC>x = (I + qQ) 
where / is the identity matrix. The solution to (|3.9|) is simply, 



N(k,x) = N(k,0) 
N(-k,x) = N(k,0) 



-q(x-L) + l 

1 + qL 
-qjx - L) 
1 + qL 



(3.17) 
(3.18) 



An example of spatial variation of the wave spectrum from x = to x = L is shown in 
Figure |21 for U = 0, and a uniform (white) incident spectrum. The reflected wave energy 



Current effects on scattering of surface gravity waves by bottom topography 17 

A(c) ' 




(b) 

reffee/ea' 
wave s/?ecfrum 




/ra/js/m/fed 
wave s/?ec/rum 



-1 J 1 

2 ^ o 

Figure 2. Bottom spectrum and evolution of a surface wave spectrum along a field of sinusoidal 
bars for U = 0, b = 0.05 m, H = 0.156 m, so that n = 6/// = 0.32, and l = 2n, n = 4, so 
that L — 4 m (bottom shown in figure 1). (a) square root of the bottom spectrum, (b) and 
(c) normalized square root wave spectrum upwave (at x < 0) and downwave (at x > L) of the 
bars, respectively. The incident spectrum (k > at x = 0) is specified to be white (unform in 
wavenumbers). 



(at k < in figure [21 a) compensates the loss of energy in the transmitted spectrum (at 
k > in figured 6). 

For k — 1/2, in the limit of small bar amplitudes, and replacing (|3.4f> in (|3.15fl yields 

Rs - («L)"° + OfaL) - 2tg f s ^ {2m) + OfaL) (3.19) 

which is identical to Mei's (1985) equation (3.21)-(3.22) for exact resonance, in the limit 
of qL <C 1, and also converges to the result of Davies & Heathershaw (1984) for that same 
limit. For large bar amplitudes, the reflection is significant if the bars occupy a length L 
longer than the localization length 1/q. However, the reflection coefficient for the wave 

1/2 

amplitude only increases with L as [qL/(l + qL)] , which is slower than the exponential 
asymptote given by Mei (1985) for sinusoidal bars, and predicted by (Belzons et al. 1988) 
from the lowest-order theory applied to a random bottom. The present inclusion of the 
correlations of second-order and zeroth order terms may be thought as the representation 
of multiple reflections that tend to increase the penetration length in the random medium. 

A deeper understanding of this question is provided by the comparison of numerical 
estimations of the reflection coefficients for the wave amplitudes R. A benchmark esti- 
mation for linear waves is provided by the step-wise model of Rey (1995) using integral 
matching conditions for the free propagating waves and three evanescent modes at the 
step boundaries. This model is known to converge to the reflection coefficents given by 
an exact solution of Laplace's equation and the boundary conditions, in the limit of an 
infinite number of steps and evanescent modes. Calculations are performed here with 70 
steps and 3 evanescent modes. These numbers are chosen because a larger number of steps 
or evanescent modes gives indistinguishable results in figureOl Results of the benchmark 
model are in good agreement with the measurements of Davies & Heathershaw (1984), 
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except for wave components for which the reflection over the beach, not included in the 
model, is comparable to the reflection over the bars. An analytical expression i?Mei is 
given by Mei (1985). R for the present second order theory is given by Rs <|3.15[1 . 

We further compare these estimates to the reflection coefficient Re,Mci that is deduced 
from the energy evolution given by Hara & Mei (1987), using the approximate solutions 
of Mei (1985, his equations 3.8-3.23). One may prefer to reformulate the energy evolution 
from the amplitude evolution equations of Kirby (1988) because he used a continuous 
water depth h = sin(mZ ), instead of Mei's h — cos(ml ) which is discontinuous at 
x = and x = L]. Yet both Mei's and Kirby's equations lead to the same energy 
exchange between the incident and reflected components. Using Mei's (1985) notations, 
the amplitudes of the incident waves, reflected waves, and bottom undulations are A = 
2<r$Q k /g, B = 2a$Q k /g, and D — — 2iG_2fc, and the 'cut-off' frequency is 

n = akD (3.20) 

2sinh(2fc£f) V ' 

The energy evolution of waves propagating over sinusoidal bars along the x-axis is given 
by Hara & Mei (1987). The reflected wave energy BB* /2 should be a solution of 

d (BB*\ c d(^) =Re(inoB * Ah (3 . 21) 



dt \ 2 J y dx\ 

where B* denotes the complex conjugate of B. This is identical to (|2.39|l for a monochro- 
matic bottom except that the imaginary part replaced by a real part. 

Equation (|3.21|) yields a corresponding energy reflection coefficient, given by the frac- 
tion of energy lost by the incoming waves, 



= -^T [ Rc{in t) B*A)dx. (3.22) 

^9 JO 



RE,Mei 

' J 3 Jo 



Simple analytical expressions can be obtained at resonance, where Mei's (1985) eq. (3.20)- 
(3.21) give, 

AB*_ = -isinh(2r(l-x/X)) 
A 2 {0) 2 cosh 2 r 

with r = ttoL/Cg, so that 

cosh2r — 1 1,2 l„i 

RemcI = — ~2 = tanh r = -i? Moi , (3.24) 

4 cosh t 2 I 

and 

-R.E, Mei = 2 _1 / 2 i?Mei- (3.25) 

It is not surprising that the energy transfer thus computed differs from the energy com- 
puted from the amplitude evolution equations. This is typical of small perturbation 
methods, and was discussed by Hasselmann (1962), among others. Yet, it is remarkable 
that the ratio of the two is exactly one half. The transfer of energy given by \S\qB* A in 
H3.21fl thus correspond to an amplitude reflection coefficient i?_E.Moi that is smaller by a 
factor 2 -1 / 2 , at resonance, compared to i?Mci (figure 3). This underprediction of the the 
reflexion of the energy by 13.24fl also has consequences for the analysis and calculation 
of wave set-up due to wave group propagation over a reflecting bottom. Indeed, the esti- 
mation of the scattering stress (|2.63|l . that contribute to the driving of long waves, was 

f Such a discontinuous bottom has a markedly different spectrum at low and high frequencies. 
The present theory, confirmed by calculations with Rey's (1995) numerical model, yield very 
different reflection coefficients for waves much shorter and much longer than the resonant waves 
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analyzed by Hara & Mei (1987) using a calculation similar to Ij3.24(l . which is a factor 
2 too small. This may explain, in part, their under-prediction of the observed elevation 
of the long wave travelling with the incident wave group. However, the present theory, 
compared to that of Hara & Mei (1987), is limited to small bar amplitudes, and fails to 
reproduce their observation of the transition from oscillatory to exponential decay in the 
spatial evolution of the wave amplitude. 

3.3. Effects of wave and bottom relative phases 
The energy exchange coefficient given by the source term always gives energy to the 
least energetic components (in the absence of currents) , and thus the energy evolution is 
monotonic. The action source term H2.39J) of order 77, that was neglected so far, may have 
any sign, and thus lead to oscillatory evolutions for the wave amplitudes, as predicted 
by Mei (1985) and observed by Hara & Mei (1987). At resonance, and for U = 0, it can 
be seen that the first-order energy product k^~ 2k m H^-39fl is equal to \AB*D/8, 

in the limit of a large number of bars. Based on Mei's (1985) approximate solution, in 
the absence of waves coming from across the bars, this quantity is purely real so that its 
imaginary part is zero and the corresponding reflection coefficient Rgi is zero. For U =/= 
this property remains as can be seen by replacing Mei's (1985) solution with Kirby's 
(1988). However, similar correlation terms were also neglected in the second order energy 
(Appendix B), so that the oscillations of the amplitude across the bar field, observed by 
Hara and Mei (1987) may occur due to terms of the same order as the scattering source 
term, including interactions of the sub-harmonic kind (Guazzelli et al. 1992). Further, 
the bottom-surface bispectrum in Si may become significant if there is a large amount of 
wave energy coming from beyond the bars. This kind of situation, e.g. due to reflection 
over a beach, was discussed by Yu & Mei (2000). 

In the absence of such a reflection, and away from resonance but for small values 
of the scattering strength parameter r = {qL) 1 / 2 = QoL/C g , the imaginary part of 
"^(tfc^o is an order {qL) 1 / 2 smaller than the real part and thus contributes a 

negligible amount to the reflection. 

3.4. Source term and deterministic results for sinusoidal bars 

For large bar amplitudes, such as n = b/H = 0.32 (figure 3. a), all theories with linearized 
bottom boundary conditions fail to capture the shift of the reflection pattern to lower 
wavenumbers. This effect was discussed by Rey (1992), and attributed to the non-linear 
nature of the dispersion relation and the rapid changes in the water depth. Reflection 
coefficients are still relatively well estimated. For these large amplitudes Mei's (1985) 
approximate solution is found to be more accurate at resonance compared to the source 
term. As expected from MAHR and proved here, i?Mci and Rs become identical as 
i] = b/H goes to zero (figure b). This fact provides a verification that the first order 
scattering term Si is different from Hara and Mei's (1987) energy transfer term, and only 
accounts for a small fraction of the reflection, a fraction that goes to zero as 77 — * 0. It is 
also found that for all bottom amplitudes, the source term expression provides a simple 
and accurate solution away from resonance. 

Nevertheless, the scattering source term cannot give an accurate description of the 
spatial variation of the wave amplitude over a deterministic bottom, as shown in figure 
0] This is related to the fact that, in MAHR, the present reflection coefficient was obtained 
from the theory of Pihl et al. (2002) after averaging over the auto-correlation scale of the 
bottom topography. The present theory can only provide an accurate description of the 
spatial evolution of the wave field over scales larger than this bottom auto-correlation 
distance. 
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Figure 3. Reflection coefficients for the wave amplitudes for U = 0, H = 0.156 m, /o = 2n, 
ra = 4. In (a) b — 0.05 so that rj — b/H = 0.32, corresponding to one of the experiments of 
Davies & Heathershaw (1984), and in (b), b = 0.01, so that n = b/H = 0.064. 



3.5. Effects of currents 

A prominent feature of solutions with current is the modification of the resonant condition 
from k = k' and I — 2k, to a' = a + IU and I = k + k' , discussed in detail by Kirby 
(1988). This shift was verified in the laboratory by Magne, Rey & Ardhuin (2005). The 
magnitude of the resonant peak is also largely enhanced for waves against the current, 
due to a general conservation of the action fluxes and the variation in the action transport 
velocity, from C g + U for the incident waves, to C — U for the reflected waves. Further, 
the modulation of the current and the surface elevation also introduce an additional 
scattering, via the M c term in the coupling coefficent (|2.47|) . Notations here assume that 
k is in the direction of the current and k' is opposite to the current. At resonance, in 
the limit 77 — s- 0, the amplitude reflection coefficient Rs given by l|3.15[) converges to the 
reflection coefficient given by Kirby (1988). Using our notations, he obtained 



-f?Kirby = 



a'(C. 9 + l01 1/2 tanh(g ^ (3 26) 



a (Cg> - U) 
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Figure 4. Spatial evolution of the incident and reflected wave amplitudes represented by trans- 
mission (T) and reflection (R) coefficients, in the near-resonant case U = 0, H = 0.156 m, 
lo = 2n m _1 , m — 10 bars, b — 5 cm, g = b/H — 0.12 and with a wave period T = 1.23 s. 
This situation corresponds to one of the experiments of Davies & Heathershaw (1984), and their 
measurements lie in the shaded area. 



with 

[0-0-1 {Cg + U) (Cg'-U)] 1 ' 2 
and Q, c = —M(k, k')b/ [4wF B (fc — fc')] . Our amplitude reflection coefficient Rs is esti- 
mated with the approximation e qC>L = (/ + qQ) L + O ((gi) 2 ), so that, to first order in 
qL, 

-a'C g C'qL^ l/2 



Rs 



G 

Replacing the analytical expression l|3.4|l in (|3.11|) yields 



(3.28) 



Rsa ^MM, 



which is clearly identical to (|3.2(il) at first order in qL. 

For finite values of qL, the reflection coefficient (|3.15fl corresponding to the solution 

of H3.9fl is obtained by calculating the proper matrix exponential. Anticipating oceano- 

graphic conditions with a water depth of 20 m, a strong 2 m s _1 current corresponds 

to a Froude number of 0.17 only. For such a low value of Fr in the context of Davies & 

Heathershaw's (1984) laboratory experiments, the convergence of the present theory and 

that of Kirby (1988) is illustrated in figure 5. The reflection coefficient is largely increased 

for following currents due to the general conservation of the wave action flux. In that case 

1/2 

R is enhanced by the factor {a(Cg + U)/ [a'(Cg' - E/)]} ■ The overall increase in R for 
following waves amounts to about 60% at Fr = 0.17, for the laboratory sinusoidal bars 
of Davies & Heathershaw (1984) shown before (figure 3), with a reflected wave energy 
multiplied by a factor 2.5, compared to the case without current. For this mild current 
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Figure 5. Amplitude reflection coefficients for monochromatic waves over sinusoidal bars for 
the same settings as in figure 3, with a following (left) or opposing (right) current of magnitude 
U — 0.2 m s _1 . For reference the reflection coefficient without current, as given by the exact 
model of Rey (1995), is also shown. The position of the resonant wavenumber is indicated with 
the grey vertical dash-dotted line. 



the contribution of the current fluctuation to the coupling coefficient is small, with a 
maximum increase of 16% on the action reflection coefficent, 8% for the wave amplitude. 
However, for larger Froude numbers, this additional scattering may become significant 
as illustrated by figure El The present theory and that of Kirby (1988) agree reasonably 
well for finite values of rj, and we thus expect the source term to represent accurately the 
scattering of waves over bottom topographies in cases of uniform currents. 

For m = 4 sinusoidal bars, the energy reflection coefficients was found to be within 
10% of the exact solution for over 90% of the wavenumber range shown in figure 3, for 
77 < 0.1 and Fr = 0, and this conclusion is expected to hold for Fr < 0.2, given the 
agreement with Kirby's (1988) approximate solution. This accuracy is twice better than 
what was found for a rectangular step with Fr — (MAHR). The present method has the 
advantage of a large economy in computing power. This method is also well adapted for 
natural sea beds, for which continuous bathymetric coverage is only available in restricted 
areas, and thus only the statistical properties of the bottom topography are accessible, 
assuming homogeneity. 



4. Scattering with current on a realistic topography 

4.1. Sandwaves in the North Sea 

A real ocean topography, at least on the continental shelf, generally presents a continuous 
and broad bottom elevation spectrum. The effects of a mean current on wave scattering 
are now examined using a bottom spectrum estimated from a detailed bathymetric survey 
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Figure 6. Amplitude reflection coefficients for monochromatic waves over sinusoidal bars for the 
same settings as in figure 3 and 4, with a stronger following current of magnitude U = 0.6 m s _1 . 
The position of the resonant wavenumber is indicated with the grey vertical dash-dotted line. The 
vertical dashed line corresponds to the wavenumber for which Cg' — U. For larger wavenumbers 
the reflected waves are blocked and cannot propagate against the current. 

of an area centered on the crest of a sand dune, in the southern North Sea (figure 0. In 
this region, tidal currents are known to generate a wide array of bedforms, from large 
scale tidal Banks to sand dunes and sand waves (e.g. Dyer & Huntley 1999; Hulscher & 
van den Brink 2001). Although sand dunes present a threat to navigation and are closely 
monitored (Idier et al. 2002), dunes are much larger than typical wind sea and swell 
wavelengths. These dunes, however, are generally covered with shorter sandwaves. In the 
surveyed area the sandwaves have a peak wavelength of 250 m, and an elevation variance 
of 1.7 m 2 , which should lead to strong oblique scattering of waves with periods of 10 s and 
longer. Over smaller areas of 3 by 3 km the variance can be as large as 3.3 m 2 with a better 
defined spectral peak, so that our chosen spectrum is expected to be representative of the 
entire region, including high and low variances on dunes crests and troughs, respectively. 
The southern North Sea is also known for the attenuation of long swells, generated in 
the Norwegian Sea. This attenuation has been generally attributed to the dissipation of 
wave energy by bottom friction (Weber 1991). 

The bottom spectrum of the chosen area, like the spectra that were obtained by AH 
from the North Carolina shelf, rolls off sharply at high wavenumbers, typically like Z~ 3 for 
the directionally-integrated bottom spectrum F B2D , and proportional to l~ A for the full 
spectrum F B . Here the maximum variance is found for bottom wavelengths of the order 
of or larger than 250 m (figure [7J. For a typical swell period of 10 s, this corresponds to 
2 times the wavelength in 20 m depth, and thus a rather small scattering angle, 30° off 
from the incident direction. Swells propagating from a distant storm, with fixed absolute 
frequency u> = a + k-U, should be reflected by bottom undulations with widely different 
variances as the current changes. 

Given this bottom spectrum and the mean water depth, simple solutions are available 
for uniform conditions, because the scattering source term is a linear function of the 
directional spectrum at a given value of the absolute frequency u (see AH for numerical 
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Figure 7. (a) high-resolution bathymetry of a sand wave field in the southern North Sea with 
depths relative to chart datum, and (b) corresponding bottom elevation spectrum with contour 
values representing log 10 (47t 2 _F s ). The locus of the interacting bottom and surface wave com- 
ponents are indicated for a 12.5 s waves from the North- East in 25 m depth, with U — (middle 
circle), U = 2 m s _1 (smaller ellipse), and U = —2 m s _1 (larger ellipse), U is positive from 
the North-East, (c) Direction-integrated bottom variance spectra from the North Carolina shelf 
and the southern North Sea. Vertical lines indicate k/l ratios and incident resonant directions 
0i, assuming an incident wave field of 12.5 s period in 25 m depth and bedforms parallel to the 
y-axis. For such bedforms, the angle between incident and scattered waves is 180° — 20 1. 

methods). We consider the wave directional spectrum for a frequency /o and discretize 
it in N a directions. This spectrum is thus a vector E in a space with N a dimensions. 
The square matrix S such that dE/df = SE is symmetric and positive, and can thus 
be diagonalized, which gives N a eigenvalues A„ and corresponding eigenvectors V„, such 
that SV„ = A„ V„. Thus the time evolution is easily obtained by a projection of E on the 
basis {V„, 1 ^ n ^ N a }, giving a decomposition of E in elementary components. Each 
of these components of the directional spectrum decays exponentially in time, except 
for the isotropic part of the spectrum which remains constant because that eigenvector 
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corresponds to A = 0. The eigenvalues thus give interesting timescales for the evolution 
of the spectrum toward this isotropic state, with a half-life time of each eigenvector given 
by - In2/A„. 

Numerical results are shown here for a mean water depth of 20 m, in order to make the 
result more visible. For that depth, waves with a period T = 10 s have a dimensionless 
depth kH = 1.04, which is close the value for which the coupling coefficient Mb is 
maximum (AH). As a result, scattering is probably stronger than in real conditions 
where the mean water depth is 30 m. The following results should still provide some 
understanding of the likely real effects, at least for larger wave periods with similar 
values of kH. Without current, if kH is kept constant, the magnitude of the coupling 
coefficient K(k,k',H) decreases like H~ 9 / 2 (AH), but it is compounded by a higher 
bottom elevation spectral density for small values of k. For back-scattering, the bottom 
wavenumbers are generally in the range where the bottom spectrum rolls off like l~ A 
(figure 0). Therefore, for these back-scattering directions, the evolution time scale of 
waves with the same value of kH, e.g. T — 11.2 s in 25 m depth or T = 13.2 s in 
35 m depth, is larger by a factor (25/20) 1/2 ~ 1.1 or (35/20) 1/2 ~ 1.3, respectively. For 
incident wave and scattering directions for which the bottom spectrum is more uniform 
and does not compensate for the reduction in the coupling coefficient, such as forward 
scattering of waves from the North- West, the time scales increase by (25/20) 9 / 2 ~ 3.7 or 
(35/20) 9 / 2 ~ 77, respectively. 

With N a = 120, corresponding to a directional resolution of 3°, figure 8 shows that 
the shortest time scales (large negative values of A n ) correspond to directional spectra 
(eigenvectors) with strong local variations. These eigenvectors are thus associated with 
scattering at small oblique angles (forward scattering). Only the last 10 eigenvalues have 
a rather broad support, corresponding to scattering at much larger angles. Besides, the 
strongest scattering corresponds to a half-life time of 430 s, and mostly affects waves from 
the North- West or South-East, i.e. propagating in a direction along the sandwave crests. 
The timescale for waves from the North-East or South- West is about five times larger (the 
corresponding range of indices is 80 < n < 110). The n — 118 eigenvector corresponds 
to an exchange of wave energy between waves travelling in opposite directions across the 
sandwaves, but the corresponding half-life is of 3 hours and 15 minutes. Similar results 
were found for N a — 180 and N a — 72 and appear little sensitive to the discretization. 

Instead of this idealized horizontally uniform situation, practical situations rather cor- 
respond to quasi-stationary conditions with spatial gradients in at least one dimension. 
In this case the simple steady solutions found above for 2D topography are not physical. 
Indeed, a 3D bottom causes scattering along the transversal direction y, and the energy 
propagating in that direction builds up slowly up to the point where it becomes as large 
as the incident wave energy. This process can take a time much longer than the typical 
duration of a storm or swell arrival, and dissipative processes are likely to be important 
as the wave energy increases (e.g. Ardhuin et al. 2003). In order to go beyond qualitative 
statements on time and spatial scales of spectral relaxation, and short of simulating an 
actual storm in two dimensions, the effects on the wave spectrum are illustrated with a 
one-dimensional model configuration. 

The source term Sbscat was introduced in the version 2.22 of the wave model WAVE- 
WATCH III (Tolman 1991, 2002), based on the wave action evolution equation l|2.57(l in 
which the time derivative on the left hand side is now a Lagrangian derivative following 
a wave packet in physical and spectral space. Bottom scattering is the only source term 
activated in the present calculation. The model was run with a spectral grid of 30 fre- 
quencies ranging from 0.04 to 0.788 Hz and a directional resolution of 3°. Unfortunately 
the model spectrum is discretized with components at fixed intrinsic frequencies a and 
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Figure 8. Eigenvalues ordered by magnitude (top), and corresponding eigenvectors (bottom 
right) of the scattering matrix S for (7 = 0, fo — 0.1 Hz, and H = 20 m. The first three and last 
three eigenvectors are shown in more detail in the bottom left. 

directions 8, which is most appropriate for other processes. Therefore a small amount of 
numerical diffusion leads to a change of action at each absolute frequencies when U ^ 0, 
and the total action is only approximately conserved in that case, with a net change of 
about 1% of the integral of the absolute value of the source term for U — ±2 m s _1 , 
and four orders of magnitudes smaller, i.e. at the round-off error level, for U — 0. We 
have chosen to show cases with significant back-scatter, corresponding to waves normally 
incident over the sandwaves. This choice also corresponds to a weaker forward scattering, 
compared to waves propagating along the the sandwave crests. 

4.2. Scattering of waves normally incident on the sandwaves 

To simplify the interpretation of the results, and the processing of the boundary con- 
ditions, a one dimensional (East- West) propagation grid is used for the computations, 
assuming that the wave field, still fully directional, is uniform in the North-South di- 
rection. The waves are propagated over a model grid 100 km long, with a mean depth 
of H = 20m, and a spatial grid step of 5 km (figure Ola). As discussed above, this wa- 
ter depth is chosen to make the result more visible, and a significant broadening of the 
incident peak with a (weaker) back-scatter of waves is also found for H — 35 m and 
f p = 0.1 Hz (not shown). A Gaussian incident surface wave spectrum is imposed, with 
a mean direction from the North-East, a narrow peak directional spread of 12°, and a 
peak frequency of 0.01 Hz (figure 0b). The source term is integrated with a time step 
of 120 s, and the advection in space uses a third order scheme with a time step of 120 s 
(Tolman 2002). 
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Figure 9. (a) Schematic of the model grid and (b) incident wave spectrum specified at point 
F. Model output is shown below for point O. Please note that waves are represented with their 
arrival direction (direction from, contrary to the standard wind sea convention). The frequency 
is the relative frequency a/2n. 



The scattering source term acts as a diffusion operator with a typical 3-lobe structure, 
negative at the peak of the wave spectrum, and positive in directions of about 30° on 
both sides of the peak. This is identical, but with a larger magnitude, to the effect 
described by AH. In general the scattering effects are relatively stronger at the lowest 
frequencies, at least in the range of frequencies used here. For still lower frequencies the 
scattering coefficient K decreases (see also AH) so that, on these spatial scales, very 
little scattering occurs for infra-gravity waves (/ < 0.05 Hz). In addition to this grazing- 
angle forward scattering, a significant back-scatter is found, in particular in the case of 
following currents. 

For an absolute wave frequency of 0.08 Hz, the curves followed by the bottom resonant 
wavenumbers are overlaid on the bottom spectrum (figuredb). The wavenumbers I along 
these curves satisfy both the relations k' + I = k and a' = a + 1 • U. Without current 
the curve is exactly a circle, and transforms to an ellipse for relatively weak currents 
(Appendix C). This approximation is used in the model to compute the source term. 
The current imposed here shifts significantly the resonant configuration for the bottom 
and surface wavenumbers. A current opposed to the waves enlarges the ellipse towards 
higher wavenumbers, while a following current will lead to a 'sampling' of shorter wave 
numbers, i.e. bottom features of larger scales. Since the bottom topography has the 
largest variance at low wavenumbers, scattering is strongest for following currents (figure 
I10|) . With our choice of parameters, there is about a factor 10 reduction in the bottom 
variance that causes backscatter as U is changed from 2 m s _1 to —2 m s _1 . Besides, the 
coupling coefficient K(k', k, H) is increased in the case of a following current, as discussed 
above for the 2D cases. 

The resulting wave spectra are also modified due to the conservation of the wave action 
flux, enhancing the reflected wave energies for U > (figure ITT|) . This effect is similar 
to what was found in the 2D cases considered above, due to the different energy flux 
velocities U + C g for the incident waves, and U — C' g for the reflected waves. In all 
cases investigated here, the narrow incident wave spectrum is significantly broadened in 
directions, and that effect is most pronounced for frequencies in the range 0.07-0.10 Hz. 
Without current or with following currents, spectra in the middle of the model domain 
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Figure 10. Computed source terms at the boundary forcing point F, (a) for U = 0, (b) for a 
following current {/ = 2 m s _1 , (c) for an opposing current U = —2 m s _1 . The frequency is the 
relative frequency / = a/2n. 



exhibit a significant level of back-scattered energy, which increases the significant wave 
height and the directional spread on the up- wave side of the sandwave field ( figure ITT|) . 
This effect should not be very sensitive to the directional spread of the incident wave 
field, because the projection of the directional spectrum on the corresponding 'smooth' 
eigenvectors of the scattering matrix (figure 8) is insensitive to local variations in the 
directional spectrum. This reflection should thus occur for a wide range of sea states. At 
the same time, the incident peak of the wave field broadens in directions as it propagates 
to the down-wave end of the model domain. This broadening is fast close the the forcing 
boundary (point F), with values of the peak frequency directional spreads ag^ p larger 
than 35° at a point 5 km inside the domain (not shown), and becomes more gradual as 
the waves propagate, due to the slower evolution of broad spectra that are associated 
with smaller eigenvalues in the scattering matrix (see also Ardhuin et al. 2003a, Ardhuin 
& Herbers 2005). It was also verified that this broadening of the main spectral peak is 
strongest for waves propagating along the main sandwave crest directions (e.g. from the 
North- West in our case) due to the larger bottom variance at I = k — k' with k ~ fc', 
resulting in a significant modification of the mean direction (Magne 2005). 

Finally, a decrease in significant wave height is found along the grid, indicating an 
attenuation due to wave-bottom scattering. In reality, bottom friction would likely induce 
a stronger decay, and that decay would be stronger than in the absence of scattering. 
Essentially the scattering increases the average time taken by wave energy to cross the 
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Figure 11. Computed wave spectra at point O, 40 km inside of the model domain, after 5 hours 
of propagation, (a) for U = 0, (b) for a following current U = 2 m s , (c) for an opposing current 
U = —2 m s _1 . The frequency is the relative frequency / = a/2n. 



domain, and, because of that longer time, bottom friction together with scattering would 
lead to a larger dissipation than friction alone (Ardhuin et al. 2003). 



5. Conclusion 

The effect of a uniform current on the scattering of random surface gravity waves 
was investigated theoretically, extending the derivations of Ardhuin & Herbers (2002). 
Wave scattering may thus be represented by a scattering source term Sbscat(fc) for each 
wave component fe, in a closed spectral action balance equation. That term gives the 
rate of exchange of wave action between wave components fc and fe' that have the same 
absolute frequency, as a result of both water depth variations on the scale of the surface 
gravity waves wavelength, and current and mean free surface inhomogeneities induced 
by the bottom topography. The exchange of action between any two wave component 
pairs k and k' is proportional to the bottom elevation spectrum at the wavenumber 
vector I = k — fe', which is characteristic of Bragg scattering. The spectral integral of the 
corresponding wave pseudo-momentum source term fc5b S cat gives a recoil force exerted 
by the bottom on the water column, in addition to the hydrostatic pressure force. 

After Magne et al. (2005a) proved that the source term was applicable to non-random 
topography and accurate in the limit of small bottom amplitudes, just like Bragg scat- 
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tering approximations for acoustic or electromagnetic waves (e.g. Elfouhaily & Guerin 
2004), it is further found here that monochromatic wave results are recovered by tak- 
ing the limit to narrow incident and reflected wave spectra. In absence of current, for 
a finite sinusoidal bottom and monochromatic waves, the reflection coefficients given by 
the source term converges to Mei's (1985) theory in the limit of the small bottom am- 
plitudes. The range of maximum reflection and the side lobe pattern of the reflection 
coefficient as a function of the incident wavenumber is thus a direct consequence of the 
shape of the bottom spectrum in that case. With this point of view, there is resonance 
at all wavenumbers but its strength is proportional to the bottom elevation variance at 
the corresponding scale. In the presence of a current, reflections converge in the same 
manner to the more general theory of Kirby (1988). In two dimensions, the main effects of 
a current is an enhancement of reflected wave amplitudes when the incident waves prop- 
agate with the current, due to a conservation of the wave action flux, and a Doppler-like 
shift of the resonant wave frequencies that undergo maximum reflection. The two scale 
approximation was found to hold very well, even for a relatively fast evolutions of the 
wave amplitudes over two wavelengths (e.g. figure 3). However, the source term does not 
give a good representation of the spatial evolution of the wave field on scales shorter 
that the bottom correlation length, nor can it give reasonable results when another wave 
train propagates from beyond the bars. In that latter lower order source term 

must be considered, and a closed action balance cannot be obtained since that extra 
term depends on the phase relationship between the incident waves, reflected waves and 
bottom undulations. 

In three dimension and over the shallow areas of the southern North Sea, where large 
sand waves are found with strong tidal currents, wave scattering is expected to be signifi- 
cant, and largely influenced by currents. Over natural topographies, the bottom typically 
de-correlates over scales shorter than the scattering-induced attenuation scales, so that 
a modification of the reflection due to a phase locking of the incident and reflected waves 
with the bottom can be neglected. The wave scattering theory presented in this paper is 
thus one more piece in the puzzle of wave propagation over shallow continental shelves, 
and this process may account for a significant part of the observed attenuation of swells in 
the southern North Sea. The representation of this phenomenon with a source term in the 
wave action balance equation is expected to be accurate in many conditions of interest. 
It is consistent with the wide use of phase-averaged models for engineering and scientific 
purposes when such large scales are involved. The alternative use of phase-resolving el- 
liptic refraction-diffraction models (e.g. Belibassakis et aZ.2001), is much more expensive 
in terms of computer resources, due to the necessity to resolve the wave phase and the 
cllipticity of the problem when back-scattering occurs. For applications to rotational cur- 
rents, the mean current U should be regarded as the wave advection velocity (Andrews 
& Mclntyre 1978, see Kirby & Chen 1989 for practical approximate expressions), but a 
detailed derivation including scattering by rotational current fluctuations should be the 
next logical extension of the present theory. This is probably achievable by coupling the 
rotational part of the flow to the irrotational part, giving a modified Bernoulli equation 
(e.g. McWilliams et al. 2004). In practice, non- homogeneities in the bottom spectrum 
will probably have to be addressed due the sharp decrease of the coupling coefficient 
with water depth, and the generally higher bottom elevation variances in the shallower 
parts of the sea floor. In particular our limited bathymetric survey shows that sandwaves 
are modulated by sand dunes, very much like short water waves are modulated by long 
waves. 
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Appendix A. Harmonic oscillator equation for the first order 
potential 

The harmonic oscillator equation i|2.32[) can be written as a linear superposition of 
equations of the type 

^A+^/i-e- 4 . (Al) 

In order to specify a unique solution to IjA If) , initial conditions must be prescribed. In 
the limit of the large propagations distances, the initial conditions contribute a negligible 
non-secular term to the solution. Following Hassclmann (1962), we choose /i(0) = and 
dfiM(O) = 0, giving, 

fl , ,x e^-e^ + ^-o/jsinM)/^ , 2 2 

fx{w,w\t) = T forw jtw , (A2) 

* , i +\ teky * smut) , 

fi(u,u';t) = -— T - ; for J = ±ui A3 
2zu/ zilj'uj 



Appendix B. Harmonic oscillator equation and energy for the second 
order potential 

Replacing <pi H2.30J) in the surface boundary condition (|2.52(l . 

(£ + **) $i > fe( * } = 9 ^ - i& ^ kH ) d ^- + 1 - VIII < ( B !) 

and conserving only the resonant terms of k , , one obtains 



dt 2 

h! ,k" 



n\Al-Ut\ 



(B2) 



with I' = (k" — k') ■ U. In order to simplify the algebra we assume that the zeroth-order 
waves are random, with no correlation between $q k and <I>q k „ unless k = ±fc" and 
s = ±s" . Thus the only contributing terms to A^o must verify k" = k. Only those terms 
are now written explicitly, the others being grouped in the '. . .'. The amplitude ^ k 
satisfies the following forced harmonic oscillator equation, 

k' 

(B3) 



^2+° 2 ) = E M2 ( fc > fc ') \°k-k' T -o--i- uy 
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This is a sum of equations of the form, 

d 2 



dt 2 



The solution f% may be written as 



<T 2 )f 2 = f 1 (a',l-U-a;t)e 11 - 



u-ut 



where 



f2,a 



te- iat ~ sin(<rf)/<7 
" 2icr p - (Z • C7 + a) 2 ] ! 



(B4) 

(B5) 
(B6) 



(2,6 



2ct' [<r' - (i • U + a)} 

e -H<r'-l-U)t i 



a 2 -(a' -I •U) 2 2a \a + (a' - I ■ U) a — (a'-l-U) 



(B7) 



The second order action contribution from correlation between the zeroth and first 
order velocity potential is given by, 



(B8) 



This correlation imposes that all non-zero terms must have k" = fe, which removes the 
'. . .' terms, so that (IB 31) becomes 



F, 



2,0, fe 

Afc 



2^M 2 (fc,fc') 



,A\ G k-k'\ ) (*t,k*o,-k) 



fe' 



Afc 



Afc 



</ 2 e iCT *)Afc, 



with 



Tit 



</*0 = ^7 V W ~ (° ~ 1 • U)] + 0(1)} • 



Taking the limit when Afc — > 0, and neglecting 0(1) terms yields 

F*(t,k)=- [ —M 2 (k,k')F B (k-k')^-^S(u;'~Lu)dk'. 
Jk' 4cr ° 



(B9) 



(BIO) 



(Bll) 



Changing the spectral coordinates from fe' to (uj',8 1 ) allows a simple removal of the 
singularity, 



*?o(t, ^ = - jT ^(k, k')F B (k - k')^f ] 



Cg' + k' ■ U/k' 



AO'. 



(B12) 



Appendix C. Resonant wavenumber configuration for U « C g 

Under the assumption U << C g , and for a current in the x direction, the resonant 
conditions 

a — a = l x U, and 
yields the following Taylor expansion to first order in a' — a, 

2" 



k' -k = (k' x -k x )-^ r + 



k 



(CI) 
(C2) 
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We define, r = k', ro — k, rcos9 — k' x , so that 

r = r + -^-(ro cos O — r cos 8), (C 3) 

and thus 

r = i — ^— s- (C4) 
1 + e cos 6 K ' 

This is the parametric equation of an ellipse of semi-major axis a, semi-minor axis 6, 

half the foci distance c, and eccentricity e, with P = ro + U/C g ro cos#o = & 2 /a, and 

e = U/Cg = c/a. The interaction between a surface wave with wavenumber k! and a 

bottom component with wavenumber I excites a surface wave with the sum wavenumber 

k = k' + I. For a fixed k and current U, in the limit of U « C g the resonant k' and / 

follow ellipses described by their polar equation (|C 4|) . that reduce to circles for U = 0. 
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